/* Solution of systems of linear equations
   Jukka Toivanen, Jan Stebel
*/

#include <math.h>
#include "solution.h"

int nprocs = 1;


int UMFPACKSolve(CSparseMat &mat, double *x, bool trans)
{
  double *null = (double *) NULL, *Ax, *sol;
//  double Control[UMFPACK_CONTROL];
  int status, n, nz, *Ap, *Ai;
  void *Symbolic, *Numeric ;
  
  n  = mat.getNrows();
  nz = mat.getNnz();
  Ap = mat.getIa();
  Ai = mat.getJa();
  Ax = mat.getV();
  
  sol = new double[n];
  
  try
  {
    status = umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    if (status < 0) throw status;
    status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null);
    if (status < 0) throw status;
    umfpack_di_free_symbolic (&Symbolic) ;
    if (trans)
    {
      status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, sol, x, Numeric, null, null);
    }
    else
    {
      status = umfpack_di_solve (UMFPACK_At, Ap, Ai, Ax, sol, x, Numeric, null, null);
    }
    umfpack_di_free_numeric (&Numeric);
    if (status < 0) throw status;
  }
  catch (int e)
  {
    //Control[UMFPACK_PRL] = 1;
    //umfpack_di_report_status(Control, e);
    delete[] sol;
    return e;
  }
  
  for (int i=0; i<n; i++)
  {
    x[i] = sol[i];
  }
  delete[] sol;
  
  return 0;
}


int solveState(CVec &r, CVec &q, double relax, int low, int ncomp) {
    CSparseMat sp;
//    cout << "compute r jacobian...\n";
    r.getJac(q(low).getInd(), ncomp, sp);
//    cout << "ok\n";
    int n = sp.getNrows();
    
    double *x = new double[n];

    for (int i=0; i<n; i++)
      x[i] = r(i).getVal();

    if (UMFPACKSolve(sp, x, false) != 0) return 1;

    for (int i=0; i<n; i++)
      q(i) = q(i)-x[i]*relax;

    delete[] x;
    return 0;
}


int solveAdjoint(CVec &r, CVar &J, CVec &q, int low, int ncomp, double *p) {
  CSparseMat sp;
  
  J.getDer(q(low).getInd(), ncomp, p);

  r.getJac(q(low).getInd(), ncomp, sp);
  
//  int res = superSolve(sp, p, true);
  
//  return res;

}


void saveMatrix(CVec &r, CVec &q, int low, int ncomp, char *file)
{
  ofstream f(file);
#ifndef AD_REV
  for (int i=0; i<r.size(); i++)
  {
    for (int j=0; j<r(i).getNder(); j++)
    {
      if (r(i).getDer(j) != 0 && r(i).getDerVar(j) >= q(low).getInd() && r(i).getDerVar(j) < q(low).getInd()+ncomp)
      {
        f << i << " " << r(i).getDerVar(j)-q(low).getInd() << " " << r(i).getDer(j) << endl;
      }
    }
  }
#endif
  f.close();
}





